Does environmental adaptation or dispersal history explain the geographical distribution of Ixodes ricinus and Ixodes persulcatus ticks in Finland?

Abstract In Finland, the distribution area of the taiga tick, Ixodes persulcatus (Schulze, 1930), is nested within a broader area of distribution of a congeneric species, the sheep tick, Ixodes ricinus (Linnaeus, 1758) (Acari: Ixodidae). We assess whether distinct environmental adaptations or dispersal history provides a more parsimonious explanation for the differences in the distributions of the two common and medically important ixodids in Finland. We used an innovative spatially constrained randomization procedure to analyze whether crowdsourced occurrence data points of the two tick species had statistically different associations with any of the 28 environmental variables. Using points of presence in a region of species co‐occurrence, we built Maxent models to examine whether environmental factors or dispersal history could explain the absence of I. persulcatus in a part of the range of I. ricinus in Finland. Five environmental variables—number of inhabitants, road length, elevation above sea level, proportion of barren bedrock and boulders, and proportion of unsorted glacial deposits—were significant at p ≤ .05, indicating greater between‐species difference in original than in the randomized data. Of these variables, only the optimum value for unsorted glacial deposits was higher for I. persulcatus than for I. ricinus. Maxent models also predicted high relative habitat suitability (suitability >80%) for I. persulcatus south of its current, sharply bounded distribution range, suggesting that the species has not fulfilled its distribution potential in Finland. The two most common and medically relevant ixodids in Finland may colonize habitats with different environmental conditions. On the contrary, the recent establishment and ongoing dispersion of I. persulcatus in Fennoscandia rather than environmental conditions cause the southernmost distribution limit of the species in Finland.


| INTRODUC TI ON
Tick-borne diseases mediated by hard ticks (Acari: Ixodidae) constitute a growing threat to public health in Europe (Marques et al., 2021). Ticks serve as vectors for a variety of pathogens (Capligina et al., 2020;Katargina et al., 2015;, of which the Borrelia burgdorferi sensu lato bacteria group (causative agents for Lyme borreliosis) and the tick-borne encephalitis virus are the most serious from a human health perspective (Bogovic & Strle, 2015). In northwestern Eurasia, the main vectors for these pathogens are two tick species: the sheep tick, Ixodes ricinus (Linnaeus, 1758), and the taiga tick, Ixodes persulcatus (Schulze, 1930) (Laaksonen et al., 2018;Marques et al., 2021).
Ixodes ricinus and I. persulcatus have been shown to co-occur in western Russia (Tokarevich et al., 2011) and in five European nations west from Russia: Estonia (Katargina et al., 2015), Finland (Laaksonen et al., 2017;Sormunen, Andersson, et al., 2020), Hungary (Hornok et al., 2020), Latvia (Capligina et al., 2020), and Sweden (Jaenson et al., 2016). Increase in abundance and expansion of the range of distribution have been reported for both species especially in the mentioned Fennoscandian regions (Finland, Sweden, and western Russia) during the past few decades (Bugmyrin et al., 2012;Tokarevich et al., 2011). Particularly the northwestern expansion of I. persulcatus that originates from western Russia appears relatively recent and rapid (Tokarevich et al., 2011), with the first established Swedish population reported from the northeastern parts of the country in 2015 (Jaenson et al., 2016). In Finland, I. persulcatus appears to have established during the late 20th century, although the lack of historical data on the occurrence of the species hampers assessing the exact timing of establishment (Öhman, 1961). For I. ricinus, historical records suggest that the species has occurred in Finland for centuries (Hulden, 2007) even if the distribution of the species was demonstrated with taxonomically and spatially reliable documentation only in the 1950s (Öhman, 1961).
The modern geographical distribution of I. persulcatus and I. ricinus in Finland has been clarified with the help of a nationwide crowdsourcing campaign (Laaksonen et al., 2017). It was found that the area of distribution of I. persulcatus is nested within the distribution area of I. ricinus ( Figure 1a). Also, the southern distribution limit of I. persulcatus was remarkably sharp and some of the patches where the species showed the highest density of occurrence lay at this distribution boundary (Figure 1b). The distribution boundary does not follow any prominent landform feature. Instead, the most obvious environmental explanation for the placement of the distribution boundary would be climate because the hemiboreal climate of southwestern Finland gradually turns to boreal climate as one advances toward north and northeast. In Eurasia, I. persulcatus has a more continental distribution than I. ricinus (Kovalev et al., 2015;

T A X O N O M Y C L A S S I F I C A T I O N
Biogeography F I G U R E 1 The occurrence points of Ixodes ricinus (a) and Ixodes persulcatus (b) in the crowdsourced data. The black polygon indicates the area of cooccurrence for the two species. Only the observations within the area of cooccurrence were used in this study.  Livanova et al., 2016). Therefore, it seems plausible that it is better adapted to dry and cold climates compared with I. ricinus (Korenberg et al., 2021). However, this climatic hypothesis is hardly the sole explanation for the southernmost distribution limit of I. persulcatus in Finland, as the species occurs and even dominates in certain regions of neighboring Estonia (Katargina et al., 2015) and Latvia (Capligina et al., 2020) where temperature sums, vegetation periods and amounts of precipitation are comparable to those in southern Finland.
In terms of biogeographical theory, our research question is an example of weighing among three alternative but not mutually exclusive explanations for spatial patterns of species distribution. Firstly, according to the niche theory (e.g., Grinnell, 1917;Hutchinson, 1957;MacArthur, 1972), every species may maintain a persistent population only in specific environmental conditions. Consequently, a viable population of a given species occurs only in such geographical areas where these environmental conditions-physical, chemical, and biological-are met (Chase & Leibold, 2003 (Hubbell, 2001).
The neutral theory states that the observed distribution of a given species results from an ongoing process of random but spatially constrained dispersal of propagules followed by a success of establishment and reproduction that is independent of the quality of the environment.
The central question in our study is whether the absence of I.
persulcatus in southern Finland is due to a difference in environmental requirements between the two tick species or whether it only shows the position of a dispersal front that eventually will advance southwards. Statistically reliable evaluation of the possible differences in environmental requirements between I. persulcatus and I.
ricinus requires that the spatial autocorrelation in their occurrence data is explicitly considered (Keitt et al., 2002;Legendre, 1993). A common approach for dealing with inflated type-I error rates due to spatial autocorrelation is to perform a restricted randomization test that preserves the spatial structure of the data. In datasets that make a regular geographical grid, this can be achieved with toroidal shift methods (Harms et al., 2001). When the sampled data are irregularly spaced, as it is in our case, the most elegant way would be the Moran spectral randomization (MRS) that generates random variables that retain the spatial structure of the original samples (Wagner & Dray, 2015). In our tick data, the MRS approach was not applicable because it cannot handle a situation in which the occurrences of two species have to be randomized simultaneously, keeping the number of co-occurrences constant. For this reason, we developed and implemented a novel spatially constrained randomization procedure to examine the environmental requirements and range of distribution of the two ixodids.
Starting from the premises laid by biogeographical theories for spatial patterns of species distribution, we examined the role of environmental factors and distribution history in explaining the distri-  (Laaksonen et al., 2018). Adults (95.5%) were the most sampled stage of tick development, followed by nymphs and larvae (4.4% and 0.1%, respectively; Laaksonen et al., 2017).
We removed individuals lacking spatial information and converted the remaining observations to presence data by removing eventual multiple individuals of the same species in one locality. In order to answer our research questions, we geographically pruned the data to the area of species co-occurrence. We delineated this area by defining a minimum bounding geometry around the loca-

| Environmental variables
Numerous environmental factors regulate the distribution of I. ricinus and I. persulcatus by influencing their host-seeking possibilities and mortality rates (Medlock et al., 2013;Sirotkin & Korenberg, 2018).
Here, we used a set of 28 environmental variables that have been shown to influence the occurrence and survival of the two ixodids.
Moreover, we included variables describing human population density and the number of summer cottages that might influence the probability of observation in crowdsourced data (Welvaert & Caley, 2016).
The explanatory variables describe soils, surface water cover, land use, climate, vegetation, human presence, and topography. A list of the explanatory data and data sources is found in Appendix.
We obtained an environmental characterization for each tick occurrence point by drawing a circular buffer with the radius of 1 km around each point and calculating statistics for each used environmental variable within the buffer (Table A1). We used the buffers to diminish the effect of probable inaccuracy in the geographical coordinates of the crowdsourced data.

| Differences in habitat preferences between I. ricinus and I. persulcatus
Our null hypothesis is that between I. ricinus and I. persulcatus there is no difference in the optima along different environmental gradients.
To test this hypothesis, we first calculated the observed difference in optimum between the species along each of the 28 environmental variables. Here, optimum was defined as the average of the values recorded for the occurrence points of the species. The spatial distributions of both I. ricinus and I. persulcatus were autocorrelated (Moran's I 0.02 and 0.01, p < .01 and < .01, respectively; tested with the package "ape," Paradis & Schliep, 2019). Accordingly, statistically reliable evaluation of the differences in environmental requirements between I. persulcatus and I. ricinus necessitates that the spatial autocorrelation in their occurrence data is explicitly considered. Moran's spectral randomization (Wagner & Dray, 2015) is a powerful method for generating randomized data while maintaining its original spatial structure. However, this randomization procedure was not applicable here because at least currently it does not allow the generation of randomized spatial structure in multi-species occurrence data that are constrained to maintain the spatial pattern of co-occurrence among the species. In our data, the observed pattern was that the two Ixodes species were never found together in the same location.
To account for the spatial autocorrelation, we developed a novel spatially constrained randomization approach that produces randomized occurrence points for the two tick species in such a way that the random points reproduce as precisely as possible the spatial pattern that we observed for I. ricinus and I. persulcatus simultaneously. A Euclidean distance matrix describes the spatial pattern of the occurrence points and, consequently, our task was to generate randomly chosen geographical coordinates of the 1232 and 1078 occurrence points so as to produce two Euclidean distance matrices that contain similar selections of distance values as possible to the two original distance matrices.
To produce such distance matrices, we randomized the tick occurrence points 10,000 times, each time separating them into two sets with 1232 and 1078 points, corresponding to N of I. ricinus and I.
persulcatus, respectively. We calculated Euclidean distance matrices based on each of the randomized sets of occurrence points (henceforth randomized distance matrices). We quantified the similarity of each randomized distance matrix to its corresponding original distance matrix by calculating an overlapping index η between the matrices (Pastore, 2018). The overlapping index is a distributioninsensitive measure of similarity between two distributions, where the similarity of distributions is defined as the proportion of overlapping area between two distributions (Pastore & Calcagni, 2019). The value of η varies from zero (no overlap) to one (complete overlap).
After each random draw, we summed the η calculated separately for I. ricinus and I. persulcatus and retained the 100 draws with highest summed η for further analysis ( Figure A1). We used the package "overlapping" (Pastore, 2018) to calculate the overlapping index.
We estimated the probability of type-I error (p) in rejecting the null hypothesis of no difference in the optimum along an environmental gradient between I. ricinus and I. persulcatus as follows. We In the hypothetical case that our 28 environmental variables were mutually uncorrelated, we would by chance expect to observe 1.4 variables with a between-species difference at a type-I risk level of 0.05 (28 × 0.05 = 1.4). Naturally, our explanatory variables are not mutually independent. To define the effective number of uncorrelated environmental dimensions, we performed a principal component analysis (PCA) on all explanatory variables used in this study, after standardizing each variable's average to zero and variance to one. Similarly, we performed a PCA with variables for which our estimated p was below 5% so that we could estimate their effective number and see whether it was greater than one would expect to observe by chance only. We used the package "factoextra" for PCA (Kassambra & Mundt, 2020).

| Absence of I. persulcatus from southern Finland
To explore whether environmental factors used in this study explain the absence of I. persulcatus from southern Finland (the second research question), we created niche models for both I.
ricinus and I. persulcatus using data from the area of species cooccurrence and predicted the probability of species occurrence in the whole of Finland. Due to the limited occurrence of either species north of the 66th parallel north (Hvidsten et al., 2014;Laaksonen et al., 2017), we used the location of the center point of Rovaniemi football stadium (approximately N° 66) as the northernmost prediction limit. The aim of this niche modeling was to explore whether the model for I. persulcatus in the area south from its distribution limit predicts a lower probability of occurrence than the model for I. ricinus in the same area. Such a difference in the model predictions would suggest that the absence of I. persulcatus in southern Finland is due to one or more of the environmental factors included in this study, whereas a lack of difference in models would suggest dispersal limitation, or an environmental factor not included in this study.
We used the maximum entropy method (Maxent) in package "dismo" (Hijmans et al., 2021) to create the niche models. Maxent is a general-purpose machine learning method, which applies the maximum entropy principle (Jaynes, 1957) for fitting a niche model in which the estimated species distribution is as close to a uniform distribution as it is possible while still being able to explain the observed dependences between occurrence points and environmental variables. Maxent is designed to be used with presence-only species observations and a set of environmental variables from the region of interest (Philips et al., 2006). Maxent predictions provide estimates of relative habitat suitability, not estimates of occupancy probability (Guillera-Arroita et al., 2014).
We fitted two kinds of models for both species: (1)  tion. We used the regularization value of 1 in each model fit. We validated the models using area under the curve value (AUC) generated using the validation data (Liu et al., 2011). The AUC values of 0.6-0.7 indicate a fair, 0.7-0.8 good, 0.8-0.9 very good, and > 0.9 excellent model classification accuracy (Swets, 1988).

| Differences in habitat preferences between I. ricinus and I. persulcatus
Out of the examined 28 environmental variables, four had an estimated p < 1% and one <5% (Figure 2, Table 1). These variables were the number of inhabitants, road length, elevation above sea level (henceforth elevation), proportion of barren bedrock and boulders (henceforth rockiness), and proportion of unsorted glacial deposits (henceforth moraine). These five variables were partially intercorrelated and a conservative estimate at 5% level of risk is that together they effectively represent between two and three independent dimensions of environmental variation (Table A2, Figure A2). This is more than one would expect to observe by chance only. A principal component analysis of the 28 environmental variables showed that effectively the environmental space can be represented by 18 independent dimensions (>95% of the total variation is captured by the first 18 principal components, Table A3). Hence, one would expect to observe an environmental difference between the two tick species along approx. One independent environmental dimension at the p ≤ .05 level (18 × 0.05 = 0.9).
Of the five environmental variables with p < .05, four showed a higher optimum value for I. ricinus than for I. persulcatus (Figure 2, Table 1). The number of inhabitants and road length can be regarded as proxies for anthropogenic factors, and they were generally higher in areas where I. ricinus were observed. However, I. persulcatus was the dominant species in major coastal cities, whereas I. ricinus prevailed in major inland cities (chi-square test of independence χ 2 = 11.3, df = 4, p = .023) ( Figure A3).

| Absence of I. persulcatus from southern Finland
All maximum entropy model predictions indicated similar patterns of relative habitat suitability for I. ricinus (Figure 3a,b) and I. persulcatus (Figure 3c,d). In predictions, relative habitat suitability for I. persulcatus was high (>0.80) also south of its current southernmost distribution limit in Finland, especially when the climate variables were excluded from the model (Figure 3d). The predictions from full models indicated low relative habitat suitability for both ixodids in lower latitudinal parts of Finland (Figure 3a,c). In both models, relative habitat suitability was higher for I. persulcatus than for I. ricinus in the coastal regions, whereas relative habitat suitability was higher for I. ricinus than for I. persulcatus in most inland regions (Figure 3e,f).

The proportion of urban area had the highest contribution in all
Maxent models, despite model type and modeled species ( In Fennoscandia, sites with a high proportion of moraine typically have low topographic positions, contrary to rocky and elevated areas (Hirvas & Nenonen, 1985). Moreover, the rocky and elevated sites often have thin or nonexistent soils and are thus poorer and more xeric than sites with the high proportion of glacial tills (Roiko-Jokela, 1980;Seibert et al., 2007). This suggests that in their area of co-occurrence, I. ricinus may colonize less fertile xeric sites than I. persulcatus.
The number of inhabitants and road length were higher in areas where I. ricinus was present than in areas where I. persulcatus Taken together, our results indicate that in the area of species co-occurrence, I. ricinus is more resistant to desiccating conditions than I. persulcatus as it inhabits sites that are more prone to drought.
Rocky and elevated sites with thin or nonexistent soils are naturally more xeric and prone to drought than low-lying sites with more finegrained soils (Hirvas & Nenonen, 1985;Roiko-Jokela, 1980;Seibert et al., 2007). Due to the urban heat island effect (UHI), urban areas typically have higher temperature than the surrounding rural areas (Kim, 1992). Increased evaporation rate caused by UHI together with effective drainage decrease moisture availability in urban areas compared with surrounding rural areas (Hage, 1975).
Our observation that I. ricinus was more often present in sites that are susceptible to drought compared with I. persulcatus may also be related to differences in temporal activity patterns between the two species. In our study region, I. ricinus has activity peaks approximately in May and August . I.
persulcatus has a single activity peak in spring/early summer, after which its activity greatly diminishes or ceases altogether ( need to conserve energy as much as possible. In xeric sites, ticks likely face additional stress due to heat and drought (Sirotkin & Korenberg, 2018). This stress may force them to search for more favorable microenvironments and/or actively uptake water, both at the expense of their energy reserves. While the same stressors influence I. ricinus, their continued activity in autumn may provide more chances to find blood meals and thus replenish energy reserves before winter dormancy, which is shorter than that of I. persulcatus.
This mechanism may lead to higher mortality rates in populations of I. persulcatus compared with I. ricinus in xeric environments.
Although the number of inhabitants and road length were higher in sites where I. ricinus occurred compared with sites where I. persulcatus was present, I. ricinus was not the dominant species in all cities.
Specifically, I. persulcatus dominated over I. ricinus in coastal cities.
The dominance of I. persulcatus over I. ricinus, has also been shown independently in the coastal city of Oulu (Pakanen et al., 2020).
We propose that the main reason as to why I. persulcatus thrives in Finnish coastal cities is due to a reverse UHI. The geographical location and spatial context of a city, such as topography and proximity of large waterbodies, may influence the warming caused by UHI (Hjort et al., 2016). Consequently, the UHI effect may reverse in highlatitude coastal cities (Suomi & Käyhkö, 2012). Moreover, the intensity of UHI correlates positively with city size (Zhou et al., 2017). Of the studied cities located on the relatively humid coastal plain of the northern Baltic Sea, Oulu is the only one with a population >70,000 (approx. 207,300 inhabitants in 2020).  (Turner & Gardner, 2015). However, variability in, e.g., microclimatic conditions is meaningful to off-host ticks that respond to environmental differences also at very fine spatial scales (Van Gestel et al., 2022).
Hence, it is possible that the occurrence of the two species is related to the environmental variables we analyzed, but the relationship was shrouded by too coarse data resolution.

| Absence of I. persulcatus from southern Finland
The Maxent models gave very similar predictions for the geographical distribution of the two species. This modeled similarity is particularly relevant in the areas where I. persulcatus was absent in our data (i.e., southern Finland). Also, the Maxent models did not indicate the existence of any belt of unsuitable habitat at the southernmost distribution limit of I. persulcatus, which could act as a dispersal barrier for the species. These observations suggest that the current south-  (Hulden, 2007;Öhman, 1961), I. persulcatus appears to have established itself relatively recently in the region (Bugmyrin et al., 2012;Jaenson et al., 2016;Laaksonen et al., 2017). Western (Tokarevich et al., 2011) and central Russia (Kovalev et al., 2015;Livanova et al., 2016) host populations of I.
persulcatus, with recent evidence of their expansion toward central and western parts of Fennoscandia (Laaksonen et al., 2017;Pakanen et al., 2020;Tokarevich et al., 2011). Both tick species rely on host animals regarding the expansion of their geographical distribution.
Birds, for example, are among the most important dispersal hosts for both species (Buczek et al., 2020;Hornok et al., 2020;Toma et al., 2021). We would expect a patchier distribution of I. persulcatus at its southernmost distribution limit should the species be vectorized mainly by animals with high mobility (Shigesada et al., 1995).
Furthermore, major southwards and westwards movements of migratory birds occur in late summer/autumn in northern Europe. As the Finnish I. persulcatus populations (and their source populations) are no longer active during this period, it is likely that they are seldom transported south or southwest by migrating birds. During the northward migration in spring, there are few places from where birds could acquire I. persulcatus-mostly inland areas in Estonia (Katargina et al., 2015) and Latvia (Capligina et al., 2020). Correspondingly, no I. persulcatus were detected in a recent survey of ticks parasitizing migrating birds in Finland . Hence, the timing and shortness of the activity period of I. persulcatus appears to contribute to limiting their distributional range in Finland, the mechanism probably being linked to the mobility of host animals. The sharp southernmost distribution limit of I. persulcatus in Finland could be due to the diffusive range expansion of the species (e.g., Shigesada et al., 1995) where ticks are transported with hosts with low mobility, such as voles (Norrdahl & Korpimäki, 1995).
On the other hand, hybridization, or competition through hybridization, i.e., converting a cumulatively increasing portion of the population of competing species into hybrids (Todesco et al., 2016) could explain the absence of I. persulcatus from southern Finland.
Hybrids formed by male I. persulcatus mating with female I. ricinus occur more commonly than vice versa (Kovalev et al., 2016).
Such higher affinity toward hybrid formation for one species may lead to the decline of the other species in areas of species cooccurrence. As southern Finland hosts abundant populations of I. ricinus (Sormunen, Andersson, et al., 2020), it appears possible that a small number of arriving I. persulcatus could be hybridized to local extinction. However, screening of approximately 3800 of the crowdsourced tick samples for I. ricinus × I. persulcatus hybrids using a duplex qPCR assay identified only 35 such hybrids (Laaksonen et al., 2018). Hence, hybridization appears to have a minor influence on our results.
The results of this study were obtained by analyzing crowdsourced data on the presence of I. ricinus and I. persulcatus in Finland. As our data were collected by citizens, the number of observations is intrinsically correlated with the number of observers (Cretois et al., 2021). This sampling bias may-in part-explain why the observations of both species were related to the occurrence of the urban area and the number of inhabitants. Together with differences in activity patterns between the two ixodids, this sampling bias may partly explain the interspecific differences in environmental associations observed in this study (Welvaert & Caley, 2016). More specifically, the longer activity period of I. ricinus compared with I. persulcatus leads to a higher probability to encounter I. ricinus than I. persulcatus. This difference may transmit to the number of observations that citizens make of each species.
The effect may be further amplified due to the concurrency of the second activity peak of I. ricinus and the main summer holiday season in Finland (Pääkkönen & Hanifi, 2011;. Despite this, the two species were observed in similar quantities from their sympatric area of occurrence (N = 1232, 1078 for I. ricinus and I. persulcatus, respectively). Still, the imbalance of species observations in our data may influence the results.
Moreover, uncertainties in, e.g., geographic information provided along the tick samples may have introduced bias in the results despite our efforts to account for the bias in the quantification of environmental variables.
We opted not to use the information on host density as it was only available for a small subset of possible hosts (over 200 reported hosts for both species; Bowmann & Nuttall, 2008;Uspensky, 2016), observations covered the study region only partially, and/or the spatial resolution of the data was too coarse. Moreover, none of the potential hosts of either species is known to have a distribution pattern that matches with the tick observations in our data. As such, we regard data on host density/occurrence as suboptimal for analyzing the niche differences of ectoparasites dependent on host animals for transportation and which may use various hosts for which occurrence data are not available. Furthermore, ecologically credible explanations for the exclusion of I. persulcatus from urban areas and rocky and elevated habitats (or vice versa for I. ricinus and moraine) based on host animal densities/presence are missing.
Should the observed differences be due to factors associated with host animals, they are more likely to be linked to the activity patterns of ticks and movement patterns of host animals during the growing season (for which data are even more sparse) than the densities of host species. Particularly, given the differences in activity patterns between the two ixodids, there is a chance that the tick species encounter different hosts and/or are transported by them differently. Thus, differences in host animal movement patterns and activity may influence the different distributions of the two ixodids in Finland.

| CON CLUS IONS
Our observations suggest that in the area of species co-occurrence,